library(AER)
library(fBasics)
library(nleqslv)
library(readstata13)
library(geosphere)

### Read in new data from Stata File "Individual Modeling"

setwd("/Users/aspearot/Dropbox/MaLT 2017/Draft/submission_draft/ReStat/resubmission 2/replication package/malt/")

sx<-read.dta13("Data/MaLT Calibration Replication.dta") 
names(sx)[1]<-"region"

## Load geolocation data to match to farmer surveys

	z_manyara<-read.csv("Data/Manyara Vil_Mkt_Turnoff_GPStoGoogle.csv",header=TRUE)
	zm<-subset(z_manyara,place=="village")
	zm$region<-"MANYARA"

	z_kili<-read.csv("Data/Kilimanjaro Vil_Mkt_Turnoff_GPStoGoogle.csv",header=TRUE)
	zk<-subset(z_kili,place=="village")

	  #remove extraneous variables in each region's file, harmonize names for variables
  	
  	zm$place<-NULL
  	zk$place<-NULL
  
  	zk$market_lat<-NULL
  	zk$market_lon<-NULL
  	zk$turnoff_lat<-NULL
  	zk$turnoff_lon<-NULL
  	zk$village_lat<-NULL
  	zk$village_lon<-NULL
  
  	zm$market_lat<-NULL
  	zm$market_lon<-NULL
  	zm$turnoff_lat<-NULL
  	zm$turnoff_lon<-NULL
  	zm$village_lat<-NULL
  	zm$village_lon<-NULL

	  names(zm)[5]<-"market"
	
	  zm$market_id<-NULL
	  zm$village_id<-NULL

	
	#merge manyara and kilimanjaro, select relevant variables, clean district and variable names
	
  	sv<-rbind(zm,zk)
	  sv<-sv[,c("region","district","ward","village_name","market","agrovet","lat","lon")]
	  sv$district<-gsub("_"," ",sv$district)
	  sv$village_name<-gsub("_"," ",sv$village_name)

###  Construct dataset for counterfactuals at the village level

	sx$unit<-paste(sx$village_name,sx$ward,sx$market,sx$district,sx$region)

	village_name<-as.character(tapply(as.character(sx$village_name),sx$unit,FUN=function(f){unique(f)}))
	farmerdata<-data.frame(village_name)
	farmerdata$ward<-as.character(tapply(as.character(sx$ward),sx$unit,FUN=function(f){unique(f)}))
	farmerdata$market<-as.character(tapply(as.character(sx$market),sx$unit,FUN=function(f){unique(f)}))
	farmerdata$region<-as.character(tapply(as.character(sx$region),sx$unit,FUN=function(f){unique(f)}))
	farmerdata$district<-as.character(tapply(as.character(sx$district),sx$unit,FUN=function(f){unique(f)}))
	
	### Sum total land, and average adoption, average expense and sum expense
	
	farmerdata$farmer_acres_land<-as.numeric(tapply(as.numeric(sx$farmer_acres_land),sx$unit,FUN=function(f){sum(f,na.rm=TRUE)}))
	farmerdata$adopt_pred<-as.numeric(tapply(as.numeric(sx$adopt_pred),sx$unit,FUN=function(f){mean(f,na.rm=TRUE)}))
	farmerdata$expense_avg<-as.numeric(tapply(as.numeric(sx$expense_pred),sx$unit,FUN=function(f){mean(f,na.rm=TRUE)}))
	farmerdata$expense_sum<-as.numeric(tapply(as.numeric(sx$expense_pred),sx$unit,FUN=function(f){sum(f,na.rm=TRUE)}))
	farmerdata$expense_pre<-as.numeric(tapply(as.numeric(sx$expense),sx$unit,FUN=function(f){mean(f,na.rm=TRUE)}))
	
	farmerdata$unit<-as.character(tapply(as.character(sx$unit),sx$unit,FUN=function(f){unique(f)}))
	farmerdata$google_vil_city_km<-as.numeric(tapply(as.numeric(sx$google_vil_city_km),sx$unit,FUN=function(f){mean(f,na.rm=TRUE)}))
	
	v<-merge(sv,farmerdata,all.x=T)
	v<-subset(v,is.na(lat)==FALSE|is.na(lon)==FALSE)	
	
	################  Construct market average of each variable  ################

	v$farmer_acres_land_market<-as.numeric(ave(as.numeric(v$farmer_acres_land),v$market,FUN=function(f){mean(f,na.rm=TRUE)}))
	v$adopt_pred_market<-as.numeric(ave(v$adopt_pred,v$market,FUN=function(f){mean(f,na.rm=TRUE)}))

	v$expense_pre_market<-as.numeric(ave(v$expense_pre,v$market,FUN=function(f){mean(f,na.rm=TRUE)}))
	v$expense_avg_market<-as.numeric(ave(v$expense_avg,v$market,FUN=function(f){mean(f,na.rm=TRUE)}))
	v$expense_sum_market<-as.numeric(ave(v$expense_sum,v$market,FUN=function(f){mean(f,na.rm=TRUE)}))	
	
	#identify which markets are sampled, and which were not

	    v$sampled<-ifelse(is.na(v$unit),0,1)

	# Load population data for both regions, and merge with farmer-level data
  	
    	popK<-read.csv("Data/Kilimanjaro VillageSize_For_Market.csv",header=TRUE)
    	popK$region<-"KILIMANJARO"
    	popM<-read.csv("Data/Manyara VillageSize_For_Market.csv",header=TRUE)	
    	popM$region<-"MANYARA"
    	
    	names(popM)[1]<-"market"
    	names(popK)[5]<-"census_households"
    	pop<-rbind(popK,popM)
	
	    v<-merge(v,pop,all.x=T)
	
	 #if population is missing, assign ward average for the population
	    
	    v$census_bothsex<-ifelse(is.na(v$census_bothsex),ave(v$census_bothsex,paste(v$region,v$ward),FUN=function(f){mean(f,na.rm=TRUE)}),v$census_bothsex)
	
##########
##########  Next, we fill in values for villages that were not in the sample and in markets that were not sampled.  This is primarily important for Manyara region.
##########
##########  To fill in values, we use averages from the closest village in the sample by straightline distance.
##########
	
  sx$id<-seq(1,nrow(sx),1)
	
	#identify observations with reported village level adoption and expenditures, and those without one of the two  
	
	    yes_v<-subset(v,is.na(adopt_pred)==FALSE&is.na(v$expense_avg)==FALSE)
	    no_v<-subset(v,is.na(adopt_pred)==TRUE|is.na(v$expense_avg)==TRUE)
	
    	no_a<-matrix(nrow=50,ncol=nrow(no_v))
    	no_f<-matrix(nrow=50,ncol=nrow(no_v))
    	no_id<-matrix(nrow=50,ncol=nrow(no_v))
    	yes_a<-matrix(nrow=50,ncol=nrow(yes_v))
    	yes_f<-matrix(nrow=50,ncol=nrow(yes_v))
    	yes_id<-matrix(nrow=50,ncol=nrow(yes_v))
	
    	#calculate distances between all villages without data and those with data
    	
    	mat<-distm(cbind(no_v$lon,no_v$lat),cbind(yes_v$lon,yes_v$lat))

    	
    for(i in 1:nrow(mat)){
    	
      #find minimum distance from village i without data to all villages with data
      
    	mindist<-which.min(mat[i,])
    	
    	#select minimum distance
    	
    	minyes_v<-yes_v[mindist,]
    
      #assign values using the minimum distance village  	
    	no_v$farmer_acres_land[i]<-minyes_v$farmer_acres_land[1]
    	no_v$adopt_pred[i]<-ifelse(is.na(no_v$adopt_pred[i]),minyes_v$adopt_pred[1],no_v$adopt_pred[i])
    	no_v$expense_avg[i]<-ifelse(is.na(no_v$expense_avg[i]),minyes_v$expense_avg[1],no_v$expense_avg[i])
    	no_v$expense_sum[i]<-ifelse(is.na(no_v$expense_sum[i]),minyes_v$expense_sum[1],no_v$expense_sum[i])
    
    	subx<-subset(sx,unit==yes_v[mindist,]$unit)
    	
    	#replicate farmer level data from the nearest village to the unsampled village.  farmers in columns, village i nrows
    	if(nrow(subx)>0){
    	  no_a[1:nrow(subx),i]<-subx$adopt_pred
    	  no_f[1:nrow(subx),i]<-subx$expense_pred
    	  no_id[1:nrow(subx),i]<-NA
    	}
    
    }
	
    #create matrix of farmers-villages for those sampled villages.  farmers in columns, village i nrows
    	
  	for(i in 1:nrow(yes_v)){
  	  
  	  subx<-subset(sx,unit==yes_v[i,]$unit)
  	  
  	    yes_a[1:nrow(subx),i]<-subx$adopt_pred
  	    yes_f[1:nrow(subx),i]<-subx$expense_pred
  	    yes_id[1:nrow(subx),i]<-subx$id
  	  
  	}	
    	
    #combine the matrices
    	
    v2<-rbind(yes_v,no_v)
    
    a2<-cbind(yes_a,no_a)  #matrix of adoption
    f2<-cbind(yes_f,no_f)  #matrix of expenditures
    id<-cbind(yes_id,no_id)  #matrix of ids

##########                         ##########
##########  Load the agrovet data  ##########
##########    Harmonize names      ##########

yK<-read.dta13("Data/Kilimanjaro agrovet_shop_level.dta")

  yK$totalfert<-yK$sellurea+yK$sellDAP+yK$sellCAN+yK$sellNPK+yK$sellOTHER
  yK$totalbio<-yK$sell_insecticide+yK$sell_herbicide+yK$sell_pesticide+yK$sell_fungicide
  yK$totalseeds<-yK$sellseeds

  yk<-yK[,c("av_ownerid","av_shopid","av_region","av_district","av_ward","av_village","av_fert_rev_16","av_fert_rev_16_USD","av_fert_rev_16_USD_w1","av_seed_rev_16","av_seed_rev_16_USD","av_seed_rev_16_USD_w1","av_fert_rev_15_USD_w1","av_seed_rev_15_USD_w1","google_vil_city_km","av_ALL_sqty_16_kg_w1","av_ALL_sqty_15_kg_w1","av_seed_sqty_16_kg_w1","av_seed_sqty_15_kg_w1","av_urea_rprice_16_USD_w1","av_urea_wprice16_USD_w1","av_urea_rprice_15_USD_w1","av_urea_wprice15_USD_w1","av_trans_cost_per50kg","av_s10","av_age","av_hh_size","totalfert","totalbio","totalseeds")]

      names(yk)[ncol(yk)-7]<-"av_urea_wprice_15_USD_w1"
      names(yk)[ncol(yk)-9]<-"av_urea_wprice_16_USD_w1"
      names(yk)[ncol(yk)-5]<-"exper"
      names(yk)[ncol(yk)-4]<-"age"
      names(yk)[ncol(yk)-3]<-"hhsize"

yM<-read.dta13("Data/Manyara agrovet_shop_level.dta")

  yM$totalfert<-yM$sellurea+yM$sellDAP+yM$sellCAN+yM$sellNPK+yM$sellOTHER
  yM$totalbio<-yM$sell_insecticide+yM$sell_herbicide+yM$sell_pesticide+yM$sell_fungicide
  yM$totalseeds<-yM$sellseeds

  ym<-yM[,c("av_ownerid","av_shopid","survey_region","av_district","av_ward","av_village","av_fert_rev_17","av_fert_rev_17_USD","av_fert_rev_17_USD_w1","av_seed_rev_17","av_seed_rev_17_USD","av_seed_rev_17_USD_w1","av_fert_rev_16_USD_w1","av_seed_rev_16_USD_w1","google_vil_city_km","av_ALL_sqty_17_kg_w1","av_ALL_sqty_16_kg_w1","av_seed_sqty_17_kg_w1","av_seed_sqty_16_kg_w1","av_urea_rprice_17_USD_w1","av_urea_wprice_17_USD_w1","av_urea_rprice_16_USD_w1","av_urea_wprice_16_USD_w1","av_trans_cost_per50kg","av_s10","av_age","av_hh_size","totalfert","totalbio","totalseeds")]

      names(ym)[ncol(ym)-5]<-"exper"
      names(ym)[ncol(ym)-4]<-"age"
      names(ym)[ncol(ym)-3]<-"hhsize"
    
    	names(ym)[13]<-"av_fert_rev_15_USD_w1"
    	names(ym)[14]<-"av_seed_rev_15_USD_w1"
    
    	names(ym)[length(names(ym))-7]<-"av_urea_wprice_15_USD_w1"
    	names(ym)[length(names(ym))-8]<-"av_urea_rprice_15_USD_w1"
    
    	names(ym)[length(names(ym))-11]<-"av_seed_sqty_15_kg_w1"
    	names(ym)[length(names(ym))-13]<-"av_ALL_sqty_15_kg_w1"
    
    	names(yk)[3:6]<-c("av_region","av_district","av_ward","av_village")
    	names(yk)<-gsub("_16","",names(yk))
    	names(yk)<-gsub("16","",names(yk))
    
    	names(ym)[3:6]<-c("av_region","av_district","av_ward","av_village")
    	names(ym)<-gsub("_17","",names(ym))
    	names(ym)<-gsub("17","",names(ym))

  # Merge files  	
    	
	  y<-rbind(yk,ym)

	# Keep only agrovet locations with positive revenues and reported retail prices for UREA
	
	  y<-subset(y,av_fert_rev>0&is.na(av_urea_rprice_USD_w1)==FALSE)
	
######## 
########  Save files for future use
######## 
	
	write.csv(v2,"Data/Temp/Farmer_Replication.csv",row.names=FALSE)	
	write.csv(y,"Data/Temp/Agrovet_Replication.csv",row.names=FALSE)
	
	write.csv(a2,"Data/Temp/Adoption_Replication.csv",row.names=FALSE)	
	write.csv(f2,"Data/Temp/Expenditure_Replication.csv",row.names=FALSE)
	write.csv(id,"Data/Temp/ID_Replication.csv",row.names=FALSE)
	
	z<-read.dta13("Data/PooledRegions ViltoVil_KMCOST.dta")
	write.csv(z,"Data/Temp/Costs_Distances_Replication.csv",row.names=FALSE)	

	z1<-read.dta13("Data/PooledRegions ViltoVil_KMCOST_share.dta")
	write.csv(z1,"Data/Temp/Costs_Shares_Replication.csv",row.names=FALSE)	
	
	######## 
	########  Load data again (you can start here if you don't want to reconstruct the basic files)
	######## 
	
	library(AER)
	library(nleqslv)
	library(readstata13)
	library(matrixStats)
	library(estimatr)
	
	v2<-read.csv("Data/Temp/Farmer_Replication.csv",header=TRUE)	
	
	y<-read.csv("Data/Temp/Agrovet_Replication.csv",header=TRUE)
	z<-read.csv("Data/Temp/Costs_Shares_Replication.csv",header=TRUE)	
	w<-read.csv("Data/Temp/Costs_Shares_Replication.csv",header=TRUE)	
	
	a2<-as.matrix(read.csv("Data/Temp/Adoption_Replication.csv",header=TRUE))	
	f2<-as.matrix(read.csv("Data/Temp/Expenditure_Replication.csv",header=TRUE))		
	id<-as.matrix(read.csv("Data/Temp/ID_Replication.csv",header=TRUE))		
	
	av_village<-as.character(tapply(as.character(v2$village_name),paste(v2$village,v2$ward,v2$district),FUN=unique))
	av_ward<-as.character(tapply(as.character(v2$ward),paste(v2$village,v2$ward,v2$district),FUN=unique))
	av_district<-as.character(tapply(as.character(v2$district),paste(v2$village,v2$ward,v2$district),FUN=unique))
	av_market<-as.character(tapply(as.character(v2$market),paste(v2$village,v2$ward,v2$district),FUN=unique))
	av_pop_market<-as.numeric(tapply(v2$census_bothsex,paste(v2$village,v2$ward,v2$district),FUN=sum,na.rm=TRUE))
	popframe<-data.frame(av_village,av_ward,av_district,av_market,av_pop_market)
	
	#Merge populations to agrovet data
  	
  	y<-merge(y,popframe,all.x=T)
  	
  	v2$census_bothsex<-v2$census_bothsex/sum(v2$census_bothsex,na.rm=TRUE)
  	
  	z$same<-ifelse(z$survey_region_O==z$survey_region_D&z$survey_district_O==z$survey_district_D&z$survey_ward_O==z$survey_ward_D&z$survey_village_O==z$survey_village_D,1,0)
	

  # Calculate rural and main-road distances
  	
  	z$DIST_km_direct_rural<-z$DIST_km_direct*z$rural_km_share	
  	z$DIST_km_direct_main<-z$DIST_km_direct*(1-z$rural_km_share)
  	
  #Ensure factors are removed from variables
  	
	v2$region<-as.character(v2$region)
	v2$district<-as.character(v2$district)
	v2$market<-as.character(v2$market)
	v2$ward<-as.character(v2$ward)
	v2$village_name<-as.character(v2$village_name)

	y$av_region<-as.character(y$av_region)
	y$av_district<-as.character(y$av_district)
	y$av_ward<-as.character(y$av_ward)
	y$av_village<-as.character(y$av_village)
	
	for(i in 1:10){
		z[,i]<-as.character(z[,i])
	}

	###  Construct bilateral matrix of distance.  J (rows) will be agrovets, and I (columns) will be villages

	J<-nrow(y)
	I<-nrow(v2)

	d<-matrix(nrow=J,ncol=I) # distance matrix
	r<-matrix(nrow=J,ncol=I) # rural share matrix

	
	#construct distance matrices
	
	for(i in 1:I){
	
		sz<-subset(z,survey_region_O==v2$region[i]&survey_district_O==v2$district[i]&survey_ward_O==v2$ward[i]&survey_village_O==v2$village[i])
	
		for(j in 1:J){	
	
			ssz<-subset(sz,survey_region_D==toupper(y$av_region[j])&survey_district_D==y$av_district[j]&survey_ward_D==y$av_ward[j]&survey_village_D==y$av_village[j])
	
			if(nrow(ssz)>0){
					d[j,i]<-ssz$DIST_km_direct[1]
					r[j,i]<-ssz$rural_km_share[1]				
				}
			if(nrow(ssz)==0){
					d[j,i]<-NA
					r[j,i]<-NA				
				}	
		}
		if(i%%100==0){print(i)}
	}
		
	minmat<-matrix(ncol=1,nrow=ncol(d))
	
	for(i in 1:ncol(d)){
	  minmat[i,1]<-min(d[,i],na.rm=TRUE)
	}
	
	rowrem<-grep(1,ifelse(rowSums(is.na(d))>500,1,0))  #agrovets with incomplete village linkages - 500 is just an arbitrarily large number chosen based on inspection of the data
	colrem<-grep(1,ifelse(colSums(is.na(d))>300,1,0))  #villages with incomplete agrovet linkages - 300 is just an arbitrarily large number chosen based on inspection of the data

	d2<-d[-rowrem,-colrem]  	#remove incomplete villages and agrovets
	r2<-r[-rowrem,-colrem]

	# farmer level adoption, expenditures, and IDs for counterfactuals
	
	a2<-a2[,-colrem]
	f2<-f2[,-colrem]
	id<-id[,-colrem]
	
	adopt_mat<-a2
	omega_mat<-f2  #assigned initially as expenditures, but make into expenditure shares
	
	for(i in 1:nrow(omega_mat)){
	  omega_mat[i,]<-omega_mat[i,]/colSums(f2,na.rm=TRUE)
	}

  #Caclulate aggregate expenditures at the village level by taking the average within the sample and multiplying by village size

	E_fert<-v2$expense_avg*v2$census_bothsex

	#Agrovet fertilizer revenues

	V_fert<-y$av_fert_rev_USD_w1
		
	#define revenues as row vectors, and remove the agrovets for which there is no distance data for their village

		V_fert<-as.matrix(V_fert[-rowrem])

	#define adoption and expenditures as column vector.  Adoption of fertilizer includes with and without seeds.

		E_fert<-t(as.matrix(E_fert[-colrem]))
		E_fert<-ifelse(is.na(E_fert),0,E_fert)

	####################            ######################## 
	#### Calculate estimated trade costs for fertilizer ####
  ####################            ########################  

	#Load data from MNL
	xtau<-read.dta13("Data/buying_fert_MNL_binned_replication_edit.dta")
	#xtau<-read.dta13("/Users/aspearot/Dropbox/MaLT 2017/Data/MaLT Project Data/Final Analysis Datasets/buying_fert_MNL_binned_032822.dta")

	#calculate matrices of main road shares and rural road shares
	m2<-d2*(1-r2)
	s2<-d2*(r2)
	
	#calculate elasticity adjusted trade costs estimates
	mtau<-ifelse(m2,NA,NA)
	mtau<-ifelse(m2>75,exp(as.numeric(xtau[9])),mtau)
	mtau<-ifelse(m2>50&m2<=75,exp(as.numeric(xtau[8])),mtau)
	mtau<-ifelse(m2>40&m2<=50,exp(as.numeric(xtau[7])),mtau)
	mtau<-ifelse(m2>30&m2<=40,exp(as.numeric(xtau[6])),mtau)
	mtau<-ifelse(m2>20&m2<=30,exp(as.numeric(xtau[5])),mtau)
	mtau<-ifelse(m2>15&m2<=20,exp(as.numeric(xtau[4])),mtau)
	mtau<-ifelse(m2>10&m2<=15,exp(as.numeric(xtau[3])),mtau)
	mtau<-ifelse(m2>5&m2<=10,exp(as.numeric(xtau[2])),mtau)
	mtau<-ifelse(m2>0&m2<=5,exp(as.numeric(xtau[1])),mtau)
	mtau<-ifelse(m2<=0,1,mtau)
	
	rtau<-ifelse(s2,NA,NA)
	rtau<-ifelse(s2>20,exp(as.numeric(xtau[14])),rtau)
	rtau<-ifelse(s2>15&s2<=20,exp(as.numeric(xtau[13])),rtau)
	rtau<-ifelse(s2>10&s2<=15,exp(as.numeric(xtau[12])),rtau)
	rtau<-ifelse(s2>5&s2<=10,exp(as.numeric(xtau[11])),rtau)
	rtau<-ifelse(s2>0&s2<=5,exp(as.numeric(xtau[10])),rtau)
	rtau<-ifelse(s2<=0,1,rtau)	

	tau<-mtau*rtau
	
	#### Calibration Step 1 - Normalize Total fertilizer expenditures and total fertilizer revenues to one.

	E_fert2<-E_fert/sum(E_fert,na.rm=TRUE)
	V_fert2<-V_fert/sum(V_fert,na.rm=TRUE)

	#Define matrix for use in the contraction mapping
	
	Ebig<-matrix(NA,nrow=nrow(d2),ncol=ncol(d2))
	
	for(j in 1:nrow(Ebig)){
		Ebig[j,]<-E_fert2
	}
  
	Dpos<-as.matrix(ifelse(c(V_fert2)==0,0,1))

	contract_start<-function(d){
		d2<-Dpos*as.matrix(d)
		dmatpre1<-d2

		dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))	
		#define lambda for agrovet j but without agrovet j in the numerator
		mshares1<-t(t((tau))/colSums(dmat1*(tau)))
		
		#predict the delta value for agrovet j
		dnew1<-(V_fert2/(rowSums(mshares1*Ebig)))
		
		#apply the normalization of deltas	
		dnew1<-dnew1/sum(dnew1)
			
		return(dnew1)
	}

		init_d0<-V_fert2
		init_d<-V_fert2

		for(i in 1:5000){
	
			new_d<-contract_start(init_d)
			if(i%%1==0){print(c(i,round(max(abs(new_d-init_d),na.rm=TRUE),digits=3)*1000000,cor(init_d,init_d0),max(init_d)))}
			if(max(abs(new_d-init_d),na.rm=TRUE)<1e-16){
				print(c(i,round(max(abs(new_d-init_d),na.rm=TRUE),digits=3)*1000000,cor(init_d,init_d0),max(init_d)))
				break}
			init_d<-new_d
		}

	share_func<-function(d){
		#d<-init_d1
		d2<-as.matrix(d)
		dmatpre1<-d2
		
		dmatpre1<-dmatpre1/sum(dmatpre1)

		dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))
		
		m1<-dmat1*(tau)
		mshares1<-t(t(m1)/colSums(m1))

		X1<-(V_fert2-rowSums(mshares1*Ebig))

		#SSR<-sum(X3^2,na.rm=TRUE)

		return(X1)		
	}
	
	#confirm solution by showing that the distance between revenues and expected expenditures at each agrovet is very small at the solution
	
	  max(abs(share_func(init_d)))
	
		dmatpre1<-as.matrix(init_d)
		dmatpre1<-dmatpre1/sum(dmatpre1)
		dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))
		
		m1<-dmat1*(tau)
		mshares1<-t(t(m1)/colSums(m1))
				
		Phi_f<-colSums(m1)
		
		v3<-v2[-colrem,]
		v3$farmers<-16

		v3$Phi_f<-as.numeric(Phi_f)
		v3$google_vil_city_km_std<-(v3$google_vil_city_km-mean(v3$google_vil_city_km,na.rm=TRUE))/sqrt(var(v3$google_vil_city_km,na.rm=TRUE))
		
		h<-read.dta13("Data/Village Distances.dta")
		
		h$market<-NULL
		
		v3$ord<-seq(1,nrow(v3),1)
		v4<-merge(v3,h,all.x=T)
		v4<-v4[order(v4$ord),]
		
		y2<-y[-rowrem,]
		
		y2$deltas<-as.numeric(init_d)	#Quality adjusted prices, "delta", from the paper
		
    #### RUN IV REGRESSIONS

    		elasreg0<-lm(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district,y2)		
    		elas_reduced<-(-as.numeric(coef(elasreg0)[2]))
    		elas_reduced
    		
    		elasreg4<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_urea_wprice_USD_w1)+log(exper)+av_district,data=y2)
    		elas_old<-(-as.numeric(coef(elasreg4)[2]))
    		elas_old
    		
    		elasreg_IV<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_urea_wprice_15_USD_w1)+log(exper)+av_district,data=subset(y2,av_pop_market>0))
    		elas_IV<-(-as.numeric(coef(elasreg_IV)[2]))
    		elas_IV
    		
    		elasreg_IV_old<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_pop_market)+log(exper)+av_district,data=subset(y2,av_pop_market>0))
    		elas_IV_old<-(-as.numeric(coef(elasreg_IV_old)[2]))
    		elas_IV_old		
    		
    		#save for future use
    		
    		write.csv(y2,"Results/IVDataset_Replication_Update.csv",row.names=FALSE,na=".")
		
	zero<-function(l){
		ifelse(is.na(l),0,l)
	}

	infna<-function(l){
		ifelse(l==Inf,NA,l)
	}
	
	deltas_before<-y2$deltas
	price_before<-y2$av_urea_rprice_USD_w1
	
	#define elasticity used for analysis
	elas<-elas_IV
		
	#calculate main counterfactual change in trade costs
	iceberg0<-tau^(-1/elas)-1
	iceberg_half<-iceberg0/2
	tau_half<-(1+iceberg_half)^(-elas)

	#calculate counterfactual change in trade costs by main and rural roads separately
	miceberg0<-mtau^(-1/elas)-1
	miceberg_half<-miceberg0/2
	riceberg0<-rtau^(-1/elas)-1
	riceberg_half<-riceberg0/2	
	tau_mainhalf<-(1+miceberg_half)^(-elas)*(1+riceberg0)^(-elas)	
	tau_ruralhalf<-(1+miceberg0)^(-elas)*(1+riceberg_half)^(-elas)		

	minmat<-matrix(ncol=6,nrow=ncol(d2))
	
	for(i in 1:ncol(d2)){
    minmat[i,1]<-which.min(d2[,i])
    minmat[i,2]<-m2[minmat[i,1],i]
    minmat[i,3]<-s2[minmat[i,1],i]
    minmat[i,4]<-iceberg0[minmat[i,1],i]    
    minmat[i,5]<-which.max(mshares1[,i])   
    minmat[i,6]<-iceberg0[minmat[i,5],i]
	}
	
	# Calculate iceberg costs to closest agrovet villages
	
	summary(minmat[,2])
	summary(minmat[,3])
	summary(minmat[,4])
	summary(minmat[,6])
	
	# Statistic on share of respondents choosing the lowest trade cost location
	sum(minmat[,4]==minmat[,6])/nrow(minmat)
	
	# Calculate average and median iceberg costs
	
	mean(colMins(iceberg0))
	median(colMins(iceberg0))
	mean(colMins(miceberg0))
	median(colMins(miceberg0))
	mean(colMins(riceberg0))
	median(colMins(riceberg0))
	
	##### Create Dataset to compare with actual choices
	
	for(i in 1:nrow(v4)){
	  
	  ag_village<-y2$av_village
	  ag_ward<-y2$av_ward
	  ag_district<-y2$av_district
	  ag_region<-y2$av_region
	  subag<-data.frame(ag_village,ag_ward,ag_district,ag_region)
	  subag$village_name<-v4$village_name[i]
	  subag$ward<-v4$ward[i]	  
	  subag$district<-v4$district[i]
	  subag$region<-v4$region[i]
	  subag$lam<-mshares1[,i]
	  subag$weighted_dist_std<-v4$weighted_dist_std[i]
	  if(i==1){bigag<-subag} 
	  if(i>1){bigag<-rbind(bigag,subag)} 
	  
	}
	
	write.csv(bigag,"Results/Predicted_Choices.csv",col.names=TRUE,row.names=FALSE)	
	
	m1_half<-dmat1*(tau_half)
	mshares1_half<-t(t(m1_half)/colSums(m1_half))
	v4$Phi_f_half<-as.numeric(colSums(m1_half))	
	
	m1_ruralhalf<-dmat1*(tau_ruralhalf)
	mshares1_ruralhalf<-t(t(m1_ruralhalf)/colSums(m1_ruralhalf))
	v4$Phi_f_ruralhalf<-as.numeric(colSums(m1_ruralhalf))
	
	m1_mainhalf<-dmat1*(tau_mainhalf)
	mshares1_mainhalf<-t(t(m1_mainhalf)/colSums(m1_mainhalf))
	v4$Phi_f_mainhalf<-as.numeric(colSums(m1_mainhalf))	
	
	sv4<-v4[,c("region","district","ward","village_name","Phi_f_mainhalf","Phi_f_ruralhalf","Phi_f_half","Phi_f","weighted_dist_std","sampled")]
	
	sx<-read.dta13("Data/MaLT Calibration Replication.dta") 
	names(sx)[1]<-"region"
	sx$id<-seq(1,nrow(sx),1)
	
	sx<-merge(sx,sv4,all.x=T)

#################################################################################
############ ITERATE THROUGH PARAMETER ESTIMATES AND COUNTERFACTUALS ############
#################################################################################	
	
	### Calculate Mark-ups

	
	  elas0<-c(elas_IV) #three elasticity estimates
	  mult0<-c(2) #three ratios between elasticity estimates and frechet parameter epsillon
  
    count<-0	
    
    ### define results matrix-size will be the number of different counterfacultuals multiplied by the number of statistics (in our case, coefficient estimates, standard errors, and pvalues)
    size<-18
    comps<-matrix(nrow=size*length(elas0)*length(mult0),ncol=11)	  
    
    for(elas in elas0){
	  for(mult in mult0){
		
    		eps<-elas*mult  #calculate epsilon
    				
    		dmatpre1<-as.matrix(init_d)
    		dmatpre1<-dmatpre1/sum(dmatpre1)  #normalize sum of quality adjusted prices to 1 (this is not required)
    
    		#Add deltas to columns
    		dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))
    		
    		m1<-dmat1*(tau)
    		mshares1<-t(t(m1)/colSums(m1))  ##calculate lambdas
    		
    		### Calculate composite shares for markup equation
    		
        		smatpre<-mshares1*Ebig
        		smat<-smatpre/rowSums(smatpre)
    		
    		zi<-colSums(omega_mat*adopt_mat,na.rm=TRUE)
    		
    		Zmat<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
    		for(j in 1:nrow(Zmat)){Zmat[j,]<-zi}
    		
    		#calculate predicted mark-up
        		y2$mshare<-elas-elas*(eps-1)/eps*rowSums(smat*mshares1*(Zmat))
        		y2$est_markup<-1/(y2$mshare)
    		
        #calculate baseline marginal cost using mark-up and prices		
    		    y2$mc<-y2$av_urea_rprice_USD_w1/(y2$est_markup+1)
    		
    		#calculate quality terms
    		    y2$Ts<-y2$deltas/y2$av_urea_rprice_USD_w1^(-elas)
    		
    		    
    		ones<-matrix(1,ncol=ncol(tau),nrow=1)		
    	
    		## calculate counterfactual changes in trade costs
    		
        		iceberg0<-tau^(-1/elas)-1
        		iceberg_half<-iceberg0/2
        		tau_half<-(1+iceberg_half)^(-elas)
        		
        		miceberg0<-mtau^(-1/elas)-1
        		miceberg_half<-miceberg0/2
        		riceberg0<-rtau^(-1/elas)-1
        		riceberg_half<-riceberg0/2	
        		
        		tau_mainhalf<-(1+miceberg_half)^(-elas)*(1+riceberg0)^(-elas)	
        		tau_ruralhalf<-(1+miceberg0)^(-elas)*(1+riceberg_half)^(-elas)		
        		
    		#define function that solves agrovet pricing equation
    		
    		pricingfunc<-function(r){
    			
    			r<-ifelse(r==0,NA,r)  ## error wrapper that is not used
    			
          #Calculate new delta's (quality adjusted prices) at new prices
    			
      			dmatpre<-y2$Ts*r^(-elas)
      			dmatpre<-ifelse(is.na(dmatpre),0,dmatpre)
    			
    			### Calculate new lambdas and the shock
    
      			dmat<-matrix(dmatpre,nrow=nrow(newtau),ncol=ncol(newtau))
      			m<-dmat*(newtau)
      			mshares<-t(t(m)/colSums(m))
    			
    			  Shock<-colSums(m)/Phi_f
    			
    			### Calculate new adoption
    			
      			adopt_mat_new<-matrix(NA,nrow=nrow(adopt_mat),ncol=ncol(adopt_mat))
      			for(k in 1:nrow(adopt_mat)){adopt_mat_new[k,]<-Shock/(Shock+(1-adopt_mat[k,])/adopt_mat[k,])}
    			
    			### Calculate new expenditures and expenditure shares
    			
      			f2_new<-matrix(NA,nrow=nrow(f2),ncol=ncol(f2))
      			for(k in 1:nrow(f2_new)){f2_new[k,]<-f2[k,]*Shock^(1/eps)*(adopt_mat_new[k,]/adopt_mat[k,])^((eps-1)/eps)}
      			  
      			omega_mat_new<-matrix(NA,nrow=nrow(omega_mat),ncol=ncol(omega_mat))
      			for(i in 1:nrow(omega_mat_new)){
      			  omega_mat_new[i,]<-f2_new[i,]/colSums(f2_new,na.rm=TRUE)
      			} 
      
    			### Calculate composite shares for markup equation
    
      			zi_new<-colSums(omega_mat_new*adopt_mat_new,na.rm=TRUE)
      			
      			Zmat_new<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
      			for(j in 1:nrow(Zmat_new)){Zmat_new[j,]<-zi_new}			
      
    			### Calculate new aggregate expenditures
    			
      			Ebig_new<-matrix(NA,nrow=nrow(Ebig),ncol=ncol(Ebig))
      			
      			for(j in 1:nrow(Ebig_new)){
      			  Ebig_new[j,]<-colMeans(f2_new,na.rm=TRUE)*v3$census_bothsex
      			}
      			
      		### Calculate mark-ups
      			
      			smatpre_new<-mshares*Ebig_new
      			smat_new<-smatpre_new/rowSums(smatpre_new)			
      
      			mshare_new<-elas-elas*(eps-1)/eps*rowSums(smat_new*mshares*(Zmat_new))
        		newmarkup<-1/(mshare_new)				
    
    			X1<-(r-newmc*(1+newmarkup))
    			X1<-ifelse(is.na(X1),0,X1)
    			return(X1)		
    			
    		}
    
    		#this function calculates outcomes
    		
    		outcomes<-function(r,cost,tn){
    			
    			r<-ifelse(r==0,NA,r)
    			
    			#Calculate new delta's (quality adjusted prices) at new prices
      			
      			dmatpre<-y2$Ts*r^(-elas)
      			dmatpre<-ifelse(is.na(dmatpre),0,dmatpre)
      			
    			### Calculate new lambdas and the shock
      			
      			dmat<-matrix(dmatpre,nrow=nrow(newtau),ncol=ncol(newtau))
      			m<-dmat*(newtau)
      			mshares<-t(t(m)/colSums(m))
      			
      			Shock<-colSums(m)/Phi_f
    			
    			### Calculate new adoption
      			
      			adopt_mat_new<-matrix(NA,nrow=nrow(adopt_mat),ncol=ncol(adopt_mat))
      			for(k in 1:nrow(adopt_mat)){adopt_mat_new[k,]<-Shock/(Shock+(1-adopt_mat[k,])/adopt_mat[k,])}
      			
    			### Calculate new expenditures
      			
      			f2_new<-matrix(NA,nrow=nrow(f2),ncol=ncol(f2))
      			for(k in 1:nrow(f2_new)){f2_new[k,]<-f2[k,]*Shock^(1/eps)*(adopt_mat_new[k,]/adopt_mat[k,])^((eps-1)/eps)}
      			
      			omega_mat_new<-matrix(NA,nrow=nrow(omega_mat),ncol=ncol(omega_mat))
      			for(i in 1:nrow(omega_mat_new)){
      			  omega_mat_new[i,]<-f2_new[i,]/colSums(f2_new,na.rm=TRUE)
      			} 
    			
    			### Calculate new revenues
      			
      			Fi<-as.matrix(colSums(f2_new,na.rm=TRUE))
      
    			### Calculate composite shares for markup equation
      			
      			zi_new<-colSums(omega_mat_new*adopt_mat_new,na.rm=TRUE)
      			
      			Zmat_new<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
      			for(j in 1:nrow(Zmat_new)){Zmat_new[j,]<-zi_new}			
      			
    			### Calculate new aggregate expenditures
      			
      			Ebig_new<-matrix(NA,nrow=nrow(Ebig),ncol=ncol(Ebig))
      			
      			for(j in 1:nrow(Ebig_new)){
      			  Ebig_new[j,]<-colMeans(f2_new,na.rm=TRUE)*v3$census_bothsex
      			}
      			
      		### Calculate mark-ups, profits
        			
      			smatpre_new<-mshares*Ebig_new
      			Rj<-rowSums(smatpre_new)			
      			smat_new<-smatpre_new/rowSums(smatpre_new)			
      					
      			mshare_new<-elas-elas*(eps-1)/eps*rowSums(smat_new*mshares*(Zmat_new))
      			marks<-1/(mshare_new)			
    			
    			  newprofits<-Rj*(marks/(marks+1))
    
    			return(list(adopt=(adopt_mat_new),expenditures=(f2_new),markup=as.numeric(marks),sales=as.numeric(Rj),profits=as.numeric(newprofits),Phihat=as.numeric(Shock)))
    
    		}
    
    		# baseline equilibrium (if this does not converge immediately, there is an error)
    		
        		init_r<-ifelse(is.na(y2$av_urea_rprice_USD_w1),0,y2$av_urea_rprice_USD_w1)
        		newmc<-y2$mc
        		newtau<-tau
        		ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
        		adopt_mat_base<-outcomes(ans$x,y2$mc,newtau)$adopt
        		exp_mat_base<-outcomes(ans$x,y2$mc,newtau)$expenditures
        		rev_base<-outcomes(ans$x,y2$mc,newtau)$sales	
        		profit_base<-outcomes(ans$x,y2$mc,newtau)$profits		
        	  markup_base<-outcomes(ans$x,y2$mc,newtau)$markup
        	  prices_base<-ans$x
        	  shock_base<-outcomes(ans$x,y2$mc,newtau)$Phihat

        # baseline counterfactual  	  
        
        		init_r<-ifelse(is.na(y2$av_urea_rprice_USD_w1),0,y2$av_urea_rprice_USD_w1)
        		newmc<-y2$mc
        		newtau<-tau_half
        		ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
        		adopt_mat_half2<-outcomes(ans$x,y2$mc,newtau)$adopt
        		exp_mat_half2<-outcomes(ans$x,y2$mc,newtau)$expenditures
        		rev_half<-outcomes(ans$x,y2$mc,newtau)$sales
        		profit_half<-outcomes(ans$x,y2$mc,newtau)$profits
        		profit_half_noentry<-profit_half
        		markup_half<-outcomes(ans$x,y2$mc,newtau)$markup
        		prices_half<-ans$x
        		shock_half<-outcomes(ans$x,y2$mc,newtau)$Phihat
    		
    		# main roads counterfactual  	  
    		
        		newtau<-tau_mainhalf
        		ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
        		adopt_mat_mainhalf2<-outcomes(ans$x,y2$mc,newtau)$adopt
        		exp_mat_mainhalf2<-outcomes(ans$x,y2$mc,newtau)$expenditures	
        		rev_mainhalf<-outcomes(ans$x,y2$mc,newtau)$sales
        		profit_mainhalf<-outcomes(ans$x,y2$mc,newtau)$profits		
        		markup_mainhalf<-outcomes(ans$x,y2$mc,newtau)$markup
        		prices_mainhalf<-ans$x
        		shock_mainhalf<-outcomes(ans$x,y2$mc,newtau)$Phihat
        		shock_mainhalf_trade<-colSums(mshares1*newtau/tau)

    		# rural roads counterfactual  
    		
        	  newtau<-tau_ruralhalf
        		ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
        		adopt_mat_ruralhalf2<-outcomes(ans$x,y2$mc,newtau)$adopt
        		exp_mat_ruralhalf2<-outcomes(ans$x,y2$mc,newtau)$expenditures				
        		rev_ruralhalf<-outcomes(ans$x,y2$mc,newtau)$sales
        		profit_ruralhalf<-outcomes(ans$x,y2$mc,newtau)$profits		
        		markup_ruralhalf<-outcomes(ans$x,y2$mc,newtau)$markup
        		prices_ruralhalf<-ans$x
        		shock_ruralhalf<-outcomes(ans$x,y2$mc,newtau)$Phihat		
    		
      	# distribution counterfactual  
        		
        		p1<-quantile(y2$av_trans_cost_per50kg,prob=0.01,na.rm=TRUE)
        		p99<-quantile(y2$av_trans_cost_per50kg,prob=0.99,na.rm=TRUE)
        		y2$av_trans_cost_per50kg_w1<-y2$av_trans_cost_per50kg
        		summary(y2$av_trans_cost_per50kg_w1)
        		y2$av_trans_cost_per50kg_w1<-ifelse(y2$av_trans_cost_per50kg_w1>p99,p99,y2$av_trans_cost_per50kg_w1)
        		y2$av_trans_cost_per50kg_w1<-ifelse(y2$av_trans_cost_per50kg_w1<p1,p1,y2$av_trans_cost_per50kg_w1)
        		summary(y2$av_trans_cost_per50kg_w1)
        		
        		trans_glm<-glm(av_trans_cost_per50kg_w1~log(google_vil_city_km),y2,family=poisson(link="log"))
        		summary(trans_glm)
        		y2$impute_tradecost<-ifelse(is.na(y2$av_trans_cost_per50kg_w1),as.numeric(predict(trans_glm,newdata=y2,type="response")),y2$av_trans_cost_per50kg_w1)/2300
    	
        		newmc<-y2$mc
        		newmc<-y2$mc-y2$impute_tradecost/2
        		newtau<-tau
        		ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))
        		adopt_halftrans<-outcomes(ans$x,newmc,newtau)$adopt		
        	  exp_halftrans<-outcomes(ans$x,newmc,newtau)$expenditures
        	  rev_halftrans<-outcomes(ans$x,newmc,newtau)$sales
        	  profit_halftrans<-outcomes(ans$x,newmc,newtau)$profits		
        		markup_halftrans<-outcomes(ans$x,newmc,newtau)$markup
        		prices_halftrans<-ans$x
        		shock_halftrans<-outcomes(ans$x,newmc,newtau)$Phihat		
    		
        # Subsidy counterfactual
        		
        		newmc<-y2$mc*0.90
        		newtau<-tau
        		ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))
        		adopt_fisp<-outcomes(ans$x,newmc,newtau)$adopt		
        		exp_fisp<-outcomes(ans$x,newmc,newtau)$expenditures
        		rev_fisp<-outcomes(ans$x,newmc,newtau)$sales
        		profit_fisp<-outcomes(ans$x,newmc,newtau)$profits	
        		markup_fisp<-outcomes(ans$x,newmc,newtau)$markup		
        		prices_fisp<-ans$x
        		shock_fisp<-outcomes(ans$x,newmc,newtau)$Phihat			

    		
        # Create dataset of farmer-level predictions from baseline, and allcounterfactual		
        		
    		N<-50
    		
    		res<-matrix(nrow=N*nrow(v4),ncol=21)
    		
    		for(i in 1:nrow(v4)){
    		
    		  res[(1+(i-1)*N):(i*N),1]<-id[,i]
    		  res[(1+(i-1)*N):(i*N),2]<-adopt_mat[,i]
    		  res[(1+(i-1)*N):(i*N),3]<-f2[,i]
    		  #res[(1+(i-1)*N):(i*N),3]<-v4$weighted_dist_std[i]
    		  #res[(1+(i-1)*N):(i*N),4]<-v4$sampled[i]	  
    		  res[(1+(i-1)*N):(i*N),4]<-adopt_mat_half2[,i]
    		  res[(1+(i-1)*N):(i*N),5]<-exp_mat_half2[,i]	  
    		  res[(1+(i-1)*N):(i*N),6]<-adopt_mat_mainhalf2[,i]
    		  res[(1+(i-1)*N):(i*N),7]<-exp_mat_mainhalf2[,i]	  
    		  res[(1+(i-1)*N):(i*N),8]<-adopt_mat_ruralhalf2[,i]
    		  res[(1+(i-1)*N):(i*N),9]<-exp_mat_ruralhalf2[,i]	  
    		  res[(1+(i-1)*N):(i*N),10]<-adopt_halftrans[,i]
    		  res[(1+(i-1)*N):(i*N),11]<-exp_halftrans[,i]
    		  res[(1+(i-1)*N):(i*N),12]<-adopt_fisp[,i]
    		  res[(1+(i-1)*N):(i*N),13]<-exp_fisp[,i]	
    		  res[(1+(i-1)*N):(i*N),14]<-i
    		  res[(1+(i-1)*N):(i*N),15]<-shock_base[i]
    		  res[(1+(i-1)*N):(i*N),16]<-shock_half[i]
    		  res[(1+(i-1)*N):(i*N),17]<-shock_mainhalf[i]
    		  res[(1+(i-1)*N):(i*N),18]<-shock_ruralhalf[i]
    		  res[(1+(i-1)*N):(i*N),19]<-shock_halftrans[i]
    		  res[(1+(i-1)*N):(i*N),20]<-shock_fisp[i]
    		  res[(1+(i-1)*N):(i*N),21]<-v4$weighted_dist_std[i]
    		}
    		
    		res<-as.data.frame(res)
    		names(res)<-c("id","adopt_base","exp_base","adopt_half","exp_half","adopt_halfmain","exp_halfmain","adopt_halfrural","exp_halfrural","adopt_halftrans","exp_halftrans","adopt_fisp","exp_fisp","village","shock_base","shock_half","shock_mainhalf","shock_ruralhalf","shock_halftrans","shock_fisp","weighted_dist_std")
    		
    		sx2<-subset(res)
    		
    		# Run regressions of baseline, and change in log outcomes, on standardized distance
    		
    		base<-summary(lm_robust(log(adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    		half<-summary(lm_robust(log(adopt_half/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		halfmain<-summary(lm_robust(log(adopt_halfmain/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		halfrural<-summary(lm_robust(log(adopt_halfrural/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		halftrans<-summary(lm_robust(log(adopt_halftrans/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		fisp<-summary(lm_robust(log(adopt_fisp/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		
    		shockreg_base<-summary(lm_robust(log(shock_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    		shockreg_half<-summary(lm_robust(log(shock_half)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		shockreg_mainhalf<-summary(lm_robust(log(shock_mainhalf)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		shockreg_ruralhalf<-summary(lm_robust(log(shock_ruralhalf)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		shockreg_halftrans<-summary(lm_robust(log(shock_halftrans)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		shockreg_fisp<-summary(lm_robust(log(shock_fisp)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		
    		baseexp<-summary(lm_robust(log(exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    		halfexp<-summary(lm_robust(log(exp_half/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		halfmainexp<-summary(lm_robust(log(exp_halfmain/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		halfruralexp<-summary(lm_robust(log(exp_halfrural/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    		halftransexp<-summary(lm_robust(log(exp_halftrans/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    		fispexp<-summary(lm_robust(log(exp_fisp/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    
    		shockreg2_base<-summary(lm_robust(I(-log(1+adopt_base*(shock_base-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    		shockreg2_half<-summary(lm_robust(I(-log(1+adopt_base*(shock_half-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		shockreg2_mainhalf<-summary(lm_robust(I(-log(1+adopt_base*(shock_mainhalf-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		shockreg2_ruralhalf<-summary(lm_robust(I(-log(1+adopt_base*(shock_ruralhalf-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		shockreg2_halftrans<-summary(lm_robust(I(-log(1+adopt_base*(shock_halftrans-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    		shockreg2_fisp<-summary(lm_robust(I(-log(1+adopt_base*(shock_fisp-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    		
    		## Assign to results file
    				
    		comps[(size*count+1):(size*count+size),1]<-elas
    		comps[(size*count+1):(size*count+size),2]<-eps	
    		
    		comps[(size*count+1),3]<-round(base$coefficients[1,1],digits=4)
    		comps[(size*count+2),3]<-round(base$coefficients[1,2],digits=4)
    		comps[(size*count+3),3]<-round(base$coefficients[1,4],digits=4)
    		comps[(size*count+4),3]<-round(half$coefficients[1,1],digits=4)
    		comps[(size*count+5),3]<-round(half$coefficients[1,2],digits=4)
    		comps[(size*count+6),3]<-round(half$coefficients[1,4],digits=4)
    		comps[(size*count+7),3]<-round(halfrural$coefficients[1,1],digits=4)
    		comps[(size*count+8),3]<-round(halfrural$coefficients[1,2],digits=4)
    		comps[(size*count+9),3]<-round(halfrural$coefficients[1,4],digits=4)
    		comps[(size*count+10),3]<-round(halfmain$coefficients[1,1],digits=4)
    		comps[(size*count+11),3]<-round(halfmain$coefficients[1,2],digits=4)
    		comps[(size*count+12),3]<-round(halfmain$coefficients[1,4],digits=4)	
    		comps[(size*count+13),3]<-round(halftrans$coefficients[1,1],digits=4)
    		comps[(size*count+14),3]<-round(halftrans$coefficients[1,2],digits=4)
    		comps[(size*count+15),3]<-round(halftrans$coefficients[1,4],digits=4)	
    		comps[(size*count+16),3]<-round(fisp$coefficients[1,1],digits=4)
    		comps[(size*count+17),3]<-round(fisp$coefficients[1,2],digits=4)
    		comps[(size*count+18),3]<-round(fisp$coefficients[1,4],digits=4)	
    
    		comps[(size*count+1),4]<-round(base$coefficients[2,1],digits=4)
    		comps[(size*count+2),4]<-round(base$coefficients[2,2],digits=4)
    		comps[(size*count+3),4]<-round(base$coefficients[2,4],digits=4)
    		comps[(size*count+4),4]<-round(half$coefficients[2,1],digits=4)
    		comps[(size*count+5),4]<-round(half$coefficients[2,2],digits=4)
    		comps[(size*count+6),4]<-round(half$coefficients[2,4],digits=4)
    		comps[(size*count+7),4]<-round(halfrural$coefficients[2,1],digits=4)
    		comps[(size*count+8),4]<-round(halfrural$coefficients[2,2],digits=4)
    		comps[(size*count+9),4]<-round(halfrural$coefficients[2,4],digits=4)
    		comps[(size*count+10),4]<-round(halfmain$coefficients[2,1],digits=4)
    		comps[(size*count+11),4]<-round(halfmain$coefficients[2,2],digits=4)
    		comps[(size*count+12),4]<-round(halfmain$coefficients[2,4],digits=4)	
    		comps[(size*count+13),4]<-round(halftrans$coefficients[2,1],digits=4)
    		comps[(size*count+14),4]<-round(halftrans$coefficients[2,2],digits=4)
    		comps[(size*count+15),4]<-round(halftrans$coefficients[2,4],digits=4)	
    		comps[(size*count+16),4]<-round(fisp$coefficients[2,1],digits=4)
    		comps[(size*count+17),4]<-round(fisp$coefficients[2,2],digits=4)
    		comps[(size*count+18),4]<-round(fisp$coefficients[2,4],digits=4)	
    		
    		comps[(size*count+1),5]<-round(baseexp$coefficients[1,1],digits=4)
    		comps[(size*count+2),5]<-round(baseexp$coefficients[1,2],digits=4)
    		comps[(size*count+3),5]<-round(baseexp$coefficients[1,4],digits=4)
    		comps[(size*count+4),5]<-round(halfexp$coefficients[1,1],digits=4)
    		comps[(size*count+5),5]<-round(halfexp$coefficients[1,2],digits=4)
    		comps[(size*count+6),5]<-round(halfexp$coefficients[1,4],digits=4)
    		comps[(size*count+7),5]<-round(halfruralexp$coefficients[1,1],digits=4)
    		comps[(size*count+8),5]<-round(halfruralexp$coefficients[1,2],digits=4)
    		comps[(size*count+9),5]<-round(halfruralexp$coefficients[1,4],digits=4)
    		comps[(size*count+10),5]<-round(halfmainexp$coefficients[1,1],digits=4)
    		comps[(size*count+11),5]<-round(halfmainexp$coefficients[1,2],digits=4)
    		comps[(size*count+12),5]<-round(halfmainexp$coefficients[1,4],digits=4)	
    		comps[(size*count+13),5]<-round(halftransexp$coefficients[1,1],digits=4)
    		comps[(size*count+14),5]<-round(halftransexp$coefficients[1,2],digits=4)
    		comps[(size*count+15),5]<-round(halftransexp$coefficients[1,4],digits=4)	
    		comps[(size*count+16),5]<-round(fispexp$coefficients[1,1],digits=4)
    		comps[(size*count+17),5]<-round(fispexp$coefficients[1,2],digits=4)
    		comps[(size*count+18),5]<-round(fispexp$coefficients[1,4],digits=4)	
    		
    		comps[(size*count+1),6]<-round(baseexp$coefficients[2,1],digits=4)
    		comps[(size*count+2),6]<-round(baseexp$coefficients[2,2],digits=4)
    		comps[(size*count+3),6]<-round(baseexp$coefficients[2,4],digits=4)
    		comps[(size*count+4),6]<-round(halfexp$coefficients[2,1],digits=4)
    		comps[(size*count+5),6]<-round(halfexp$coefficients[2,2],digits=4)
    		comps[(size*count+6),6]<-round(halfexp$coefficients[2,4],digits=4)
    		comps[(size*count+7),6]<-round(halfruralexp$coefficients[2,1],digits=4)
    		comps[(size*count+8),6]<-round(halfruralexp$coefficients[2,2],digits=4)
    		comps[(size*count+9),6]<-round(halfruralexp$coefficients[2,4],digits=4)
    		comps[(size*count+10),6]<-round(halfmainexp$coefficients[2,1],digits=4)
    		comps[(size*count+11),6]<-round(halfmainexp$coefficients[2,2],digits=4)
    		comps[(size*count+12),6]<-round(halfmainexp$coefficients[2,4],digits=4)	
    		comps[(size*count+13),6]<-round(halftransexp$coefficients[2,1],digits=4)
    		comps[(size*count+14),6]<-round(halftransexp$coefficients[2,2],digits=4)
    		comps[(size*count+15),6]<-round(halftransexp$coefficients[2,4],digits=4)	
    		comps[(size*count+16),6]<-round(fispexp$coefficients[2,1],digits=4)
    		comps[(size*count+17),6]<-round(fispexp$coefficients[2,2],digits=4)
    		comps[(size*count+18),6]<-round(fispexp$coefficients[2,4],digits=4)	
    
    		comps[(size*count+1),7]<-round(shockreg_base$coefficients[1,1],digits=4)
    		comps[(size*count+2),7]<-round(shockreg_base$coefficients[1,2],digits=4)
    		comps[(size*count+3),7]<-round(shockreg_base$coefficients[1,4],digits=4)
    		comps[(size*count+4),7]<-round(shockreg_half$coefficients[1,1],digits=4)
    		comps[(size*count+5),7]<-round(shockreg_half$coefficients[1,2],digits=4)
    		comps[(size*count+6),7]<-round(shockreg_half$coefficients[1,4],digits=4)
    		comps[(size*count+7),7]<-round(shockreg_ruralhalf$coefficients[1,1],digits=4)
    		comps[(size*count+8),7]<-round(shockreg_ruralhalf$coefficients[1,2],digits=4)
    		comps[(size*count+9),7]<-round(shockreg_ruralhalf$coefficients[1,4],digits=4)
    		comps[(size*count+10),7]<-round(shockreg_mainhalf$coefficients[1,1],digits=4)
    		comps[(size*count+11),7]<-round(shockreg_mainhalf$coefficients[1,2],digits=4)
    		comps[(size*count+12),7]<-round(shockreg_mainhalf$coefficients[1,4],digits=4)	
    		comps[(size*count+13),7]<-round(shockreg_halftrans$coefficients[1,1],digits=4)
    		comps[(size*count+14),7]<-round(shockreg_halftrans$coefficients[1,2],digits=4)
    		comps[(size*count+15),7]<-round(shockreg_halftrans$coefficients[1,4],digits=4)	
    		comps[(size*count+16),7]<-round(shockreg_fisp$coefficients[1,1],digits=4)
    		comps[(size*count+17),7]<-round(shockreg_fisp$coefficients[1,2],digits=4)
    		comps[(size*count+18),7]<-round(shockreg_fisp$coefficients[1,4],digits=4)	
    		
    		comps[(size*count+1),8]<-round(shockreg_base$coefficients[2,1],digits=4)
    		comps[(size*count+2),8]<-round(shockreg_base$coefficients[2,2],digits=4)
    		comps[(size*count+3),8]<-round(shockreg_base$coefficients[2,4],digits=4)
    		comps[(size*count+4),8]<-round(shockreg_half$coefficients[2,1],digits=4)
    		comps[(size*count+5),8]<-round(shockreg_half$coefficients[2,2],digits=4)
    		comps[(size*count+6),8]<-round(shockreg_half$coefficients[2,4],digits=4)
    		comps[(size*count+7),8]<-round(shockreg_ruralhalf$coefficients[2,1],digits=4)
    		comps[(size*count+8),8]<-round(shockreg_ruralhalf$coefficients[2,2],digits=4)
    		comps[(size*count+9),8]<-round(shockreg_ruralhalf$coefficients[2,4],digits=4)
    		comps[(size*count+10),8]<-round(shockreg_mainhalf$coefficients[2,1],digits=4)
    		comps[(size*count+11),8]<-round(shockreg_mainhalf$coefficients[2,2],digits=4)
    		comps[(size*count+12),8]<-round(shockreg_mainhalf$coefficients[2,4],digits=4)	
    		comps[(size*count+13),8]<-round(shockreg_halftrans$coefficients[2,1],digits=4)
    		comps[(size*count+14),8]<-round(shockreg_halftrans$coefficients[2,2],digits=4)
    		comps[(size*count+15),8]<-round(shockreg_halftrans$coefficients[2,4],digits=4)	
    		comps[(size*count+16),8]<-round(shockreg_fisp$coefficients[2,1],digits=4)
    		comps[(size*count+17),8]<-round(shockreg_fisp$coefficients[2,2],digits=4)
    		comps[(size*count+18),8]<-round(shockreg_fisp$coefficients[2,4],digits=4)			
    			
    		comps[(size*count+1),9]<-round(shockreg2_base$coefficients[1,1],digits=4)
    		comps[(size*count+2),9]<-round(shockreg2_base$coefficients[1,2],digits=4)
    		comps[(size*count+3),9]<-round(shockreg2_base$coefficients[1,4],digits=4)
    		comps[(size*count+4),9]<-round(shockreg2_half$coefficients[1,1],digits=4)
    		comps[(size*count+5),9]<-round(shockreg2_half$coefficients[1,2],digits=4)
    		comps[(size*count+6),9]<-round(shockreg2_half$coefficients[1,4],digits=4)
    		comps[(size*count+7),9]<-round(shockreg2_ruralhalf$coefficients[1,1],digits=4)
    		comps[(size*count+8),9]<-round(shockreg2_ruralhalf$coefficients[1,2],digits=4)
    		comps[(size*count+9),9]<-round(shockreg2_ruralhalf$coefficients[1,4],digits=4)
    		comps[(size*count+10),9]<-round(shockreg2_mainhalf$coefficients[1,1],digits=4)
    		comps[(size*count+11),9]<-round(shockreg2_mainhalf$coefficients[1,2],digits=4)
    		comps[(size*count+12),9]<-round(shockreg2_mainhalf$coefficients[1,4],digits=4)	
    		comps[(size*count+13),9]<-round(shockreg2_halftrans$coefficients[1,1],digits=4)
    		comps[(size*count+14),9]<-round(shockreg2_halftrans$coefficients[1,2],digits=4)
    		comps[(size*count+15),9]<-round(shockreg2_halftrans$coefficients[1,4],digits=4)	
    		comps[(size*count+16),9]<-round(shockreg2_fisp$coefficients[1,1],digits=4)
    		comps[(size*count+17),9]<-round(shockreg2_fisp$coefficients[1,2],digits=4)
    		comps[(size*count+18),9]<-round(shockreg2_fisp$coefficients[1,4],digits=4)	
    		
    		comps[(size*count+1),10]<-round(shockreg2_base$coefficients[2,1],digits=4)
    		comps[(size*count+2),10]<-round(shockreg2_base$coefficients[2,2],digits=4)
    		comps[(size*count+3),10]<-round(shockreg2_base$coefficients[2,4],digits=4)
    		comps[(size*count+4),10]<-round(shockreg2_half$coefficients[2,1],digits=4)
    		comps[(size*count+5),10]<-round(shockreg2_half$coefficients[2,2],digits=4)
    		comps[(size*count+6),10]<-round(shockreg2_half$coefficients[2,4],digits=4)
    		comps[(size*count+7),10]<-round(shockreg2_ruralhalf$coefficients[2,1],digits=4)
    		comps[(size*count+8),10]<-round(shockreg2_ruralhalf$coefficients[2,2],digits=4)
    		comps[(size*count+9),10]<-round(shockreg2_ruralhalf$coefficients[2,4],digits=4)
    		comps[(size*count+10),10]<-round(shockreg2_mainhalf$coefficients[2,1],digits=4)
    		comps[(size*count+11),10]<-round(shockreg2_mainhalf$coefficients[2,2],digits=4)
    		comps[(size*count+12),10]<-round(shockreg2_mainhalf$coefficients[2,4],digits=4)	
    		comps[(size*count+13),10]<-round(shockreg2_halftrans$coefficients[2,1],digits=4)
    		comps[(size*count+14),10]<-round(shockreg2_halftrans$coefficients[2,2],digits=4)
    		comps[(size*count+15),10]<-round(shockreg2_halftrans$coefficients[2,4],digits=4)	
    		comps[(size*count+16),10]<-round(shockreg2_fisp$coefficients[2,1],digits=4)
    		comps[(size*count+17),10]<-round(shockreg2_fisp$coefficients[2,2],digits=4)
    		comps[(size*count+18),10]<-round(shockreg2_fisp$coefficients[2,4],digits=4)			
    		
    		
    		comps[(size*count+1),11]<-"base"
    		comps[(size*count+2),11]<-"base"
    		comps[(size*count+3),11]<-"base"
    		comps[(size*count+4),11]<-"half"
    		comps[(size*count+5),11]<-"half"
    		comps[(size*count+6),11]<-"half"
    		comps[(size*count+7),11]<-"rural"
    		comps[(size*count+8),11]<-"rural"	
    		comps[(size*count+9),11]<-"rural"
    		comps[(size*count+10),11]<-"main"
    		comps[(size*count+11),11]<-"main"	
    		comps[(size*count+12),11]<-"main"	
    		comps[(size*count+13),11]<-"dist"
    		comps[(size*count+14),11]<-"dist"	
    		comps[(size*count+15),11]<-"dist"
    		comps[(size*count+16),11]<-"fisp"
    		comps[(size*count+17),11]<-"fisp"	
    		comps[(size*count+18),11]<-"fisp"	
    		
    		count<-count+1
		
	  }}
				
    ### write results
       
    write.csv(comps,"Results/Replication_Log_Update.csv",row.names=FALSE)
    

    
    
    
    
    ### Welfare Calculations
    
    
    Farmer_Pre<-sum(exp_mat_base/adopt_mat_base,na.rm=TRUE)/sum(exp_mat_base/adopt_mat_base,na.rm=TRUE)/mult
    Farmer_Half<-(sum(exp_mat_half2/adopt_mat_half2,na.rm=TRUE)/sum(exp_mat_base/adopt_mat_base,na.rm=TRUE))/mult
    Farmer_MainHalf<-(sum(exp_mat_mainhalf2/adopt_mat_mainhalf2,na.rm=TRUE)/sum(exp_mat_base/adopt_mat_base,na.rm=TRUE))/mult
    Farmer_RuralHalf<-(sum(exp_mat_ruralhalf2/adopt_mat_ruralhalf2,na.rm=TRUE)/sum(exp_mat_base/adopt_mat_base,na.rm=TRUE))/mult    
    Farmer_HalfTrans<-(sum(exp_halftrans/adopt_halftrans,na.rm=TRUE)/sum(exp_mat_base/adopt_mat_base,na.rm=TRUE))/mult    
    Farmer_Fisp<-(sum(exp_fisp/adopt_fisp,na.rm=TRUE)/sum(exp_mat_base/adopt_mat_base,na.rm=TRUE))/mult        
    
    Agro_Pre<-sum(profit_base,na.rm=TRUE)/sum(rev_base,na.rm=TRUE)
    Agro_Half<-sum(profit_half,na.rm=TRUE)/sum(rev_base,na.rm=TRUE)
    Agro_MainHalf<-sum(profit_mainhalf,na.rm=TRUE)/sum(rev_base,na.rm=TRUE)
    Agro_RuralHalf<-sum(profit_ruralhalf,na.rm=TRUE)/sum(rev_base,na.rm=TRUE)
    Agro_HalfTrans<-sum(profit_halftrans,na.rm=TRUE)/sum(rev_base,na.rm=TRUE)
    Agro_Fisp<-sum(profit_fisp,na.rm=TRUE)/sum(rev_base,na.rm=TRUE)
    
    Surplus_Pre<-Farmer_Pre+Agro_Pre
    Surplus_Half<-Farmer_Half+Agro_Half
    Surplus_MainHalf<-Farmer_MainHalf+Agro_MainHalf
    Surplus_RuralHalf<-Farmer_RuralHalf+Agro_RuralHalf
    Surplus_HalfTrans<-Farmer_HalfTrans+Agro_HalfTrans
    Surplus_Fisp<-Farmer_Fisp+Agro_Fisp
    
    Farmer_Share_Pre<-Farmer_Pre/Surplus_Pre
    Farmer_Share_Half<-Farmer_Half/Surplus_Half    
    Farmer_Share_MainHalf<-Farmer_MainHalf/Surplus_MainHalf    
    Farmer_Share_RuralHalf<-Farmer_RuralHalf/Surplus_RuralHalf    
    Farmer_Share_HalfTrans<-Farmer_HalfTrans/Surplus_HalfTrans    
    Farmer_Share_Fisp<-Farmer_Fisp/Surplus_Fisp    
    
    
    
    Surplus_Half/Surplus_Pre
    Surplus_MainHalf/Surplus_Pre
    Surplus_RuralHalf/Surplus_Pre
    Surplus_HalfTrans/Surplus_Pre
    Surplus_Fisp/Surplus_Pre
    
    Farmer_Share_Pre
    Farmer_Share_Half
    Farmer_Share_MainHalf
    Farmer_Share_RuralHalf  
    Farmer_Share_HalfTrans   
    Farmer_Share_Fisp   
    
    
    
